/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s

/* Copyright (C) 2010 LinBox
 * Written by <brice.boyer@imag.fr>
 * Adapted from nullspaceMP in Storjohann and Chen's IML.
 *
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/** \file linbox/algorithms/NullspaceZ.inl
 * @brief The right or left nullspace of an integer matrix A
 * @details Provides :
 *	- the nullspace of a matrix \p A
 *	- (soon) a random vector within the nullspace of \p A
 *	.
 * @todo template by <class _IntegerRing>
 */

#ifndef __LINBOX_integer_nullspace_INL
#define __LINBOX_integer_nullspace_INL

#include "linbox/matrix/factorized-matrix.h"
#include "linbox/algorithms/matrix-hom.h"
#include "linbox/algorithms/rational-solver.h"
#include "linbox/blackbox/blas-blackbox.h"
#include "linbox/randiter/random-prime.h"
#include "linbox/integer.h"
#include "linbox/util/debug.h"
#include "linbox/algorithms/integer-tools.h"

// CRA
#include "linbox/algorithms/cra-domain.h"
#include "linbox/algorithms/cra-full-multip-fixed.h"

#include <stdlib.h>

#define _LB_MAX_NS_TRY 5
// #define __LINBOX_CUTS_HAIR_IN_FOUR

#ifdef __LINBOX_HAVE_IML
#include "iml.h"
#endif
#define __LINBOX_USE_IML
#define _LB_LOG2 0.6931471807

namespace LinBox
{ /*  useful tools */
	/// ça existe ?
	template<class Matrix>
	struct StackMat {
		Matrix & first ;
		Matrix & second;
		StackMat(Matrix & A, Matrix & B) :
			first(A),second(B)
		{}
	} ;

	/*! Extraction of submatrices according to permuations.
	 * Let \c P and \c Q be two permutations and \c A be a dense matrix.
	 * Let \c QAP be the resulting permuted matrix. For a given \c r,
	 * we want to extract the left upper \c r square matrix \c C in \c QAP and
	 * the rectangular matrix \c B formed by the remaining elements in the first \c r
	 * rows.
	 * We use the formula \f$\left(QAP\right)_{i,j} = A_{Q^{-1}(i),P(j)}\f$.
	 * @param[out] QAP such that \c QAP=pair(B,C).
	 * @param A input matrix
	 * @param Qt row permutation
	 * @param P column permutation
	 */
	template <class DenseMatrix, class DenseMatrix2, class Permutation>
	StackMat<DenseMatrix> & extractMatrix(  StackMat<DenseMatrix> & QAP,
						const DenseMatrix2 & A,
						const Permutation & Qt, const Permutation & P)
	{
		// std::cout                  << "Extracting matrix..." << std::endl;
		// A.write(std::cout          << "A is : "              << std::endl);
		// P.write(std::cout          << "P^{-1} : " , false )          << std::endl;
		// Qt.write(std::cout          << "Q^{-1} : " , false )          << std::endl;
		DenseMatrix & D = QAP.first ;
		DenseMatrix & G = QAP.second ;
		// Permutation P(Pt) ;
		// P.write(std::cout<< 'a')<<std::endl;
		// P.Invert();
		// P.Compress();
		// P.write(std::cout<< 'b')<<std::endl;
		// Pt.write(std::cout<< 'a')<<std::endl;
		size_t r = D.rowdim();
		if (!P.isIdentity()) {
			// std::cout << "P != ID" << std::endl;
			if (!Qt.isIdentity()) {
				// std::cout << "Qt != ID" << std::endl;
				for (size_t i = 0 ; i < r ; ++i){
					size_t i_r = Qt[i] ;
					size_t j = 0 ;
					for ( ; j < r ; ++j){
						D.setEntry(i,j,A.getEntry(i_r,P[j]));
					}
					for ( ; j < A.coldim() ; ++j){
						G.setEntry(i,j-r,A.getEntry(i_r,P[j]));
					}

				}
			}
			else { // Qt ==id
				// std::cout << "Qt == ID" << std::endl;
				for (size_t i = 0 ; i < r ; ++i){
					size_t j = 0 ;
					for ( ; j < r ; ++j){
						D.setEntry(i,j,A.getEntry(i,P[j]));
					}
					for ( ; j < A.coldim() ; ++j){
						G.setEntry(i,j-r,A.getEntry(i,P[j]));
					}

				}
			}
		}
		else { // P == ID
			// std::cout << "P == ID" << std::endl;
			if (!Qt.isIdentity()) {
				// std::cout << "Qt != ID" << std::endl;
				for (size_t i = 0 ; i < r ; ++i){
					size_t i_r = Qt[i] ;
					size_t j = 0 ;
					for ( ; j < r ; ++j){
						D.setEntry(i,j,A.getEntry(i_r,j));
					}
					for ( ; j < A.coldim() ; ++j){
						G.setEntry(i,j-r,A.getEntry(i_r,j));
					}

				}
			}
			else {
				// std::cout << "Qt == ID" << std::endl;
				std::cerr << "you should not be here" << std::endl;
				for (size_t i = 0 ; i < r ; ++i){
					size_t j = 0 ;
					for ( ; j < r ; ++j){
						D.setEntry(i,j,A.getEntry(i,j));
					}
					for ( ; j < A.coldim() ; ++j){
						G.setEntry(i,j-r,A.getEntry(i,j));
					}

				}
			}
		}

		// std::cout                  << "... matrix extracted." << std::endl;
		// QAP.first.write(std::cout  << "C is : "              << std::endl,false) << ';' << std::endl;
		// QAP.second.write(std::cout << "B is : "              << std::endl,false) << ';' << std::endl;

		return QAP ;
	}


#if 0
	struct ReduceForm {
		BlasMatrix<Integer>      *_C_ ;
		BlasMatrix<Integer>      *_B_ ;
		BlasPermutation<size_t>  &_P_ ;
		BlasPermutation<size_t>  &_Q_ ;
		size_t                    _r_ ;
		size_t                    _s_ ;
		ReduceForm( const BlasMatrix<Integer > & A)
		{
			typedef Modular<double> Field ;
			typedef Field::Element Element;
			size_t n = A.coldim();
			size_t m = A.rowdim();

			_P_=BlasPermutation<size_t>(n);
			_Q_=BlasPermutation<size_t>(m);

			RandomPrimeIter PrimeGen(19) ;
			Element p = (Element) PrimeGen.random_between(15);
			Field F(p);
			BlasMatrix<Element> Ap(A,F);
			LQUPMatrix<Field> LUp(F,Ap,_P_,_Q_);
			_r_ = LUp.getrank();
			_s_ = n-_r_ ;
			if (_r_ && _s_) { /* extraction */
				_P_.Compress();
				_Q_.Compress();
				if (!_Q_.isIdentity()) {
					_Q_.Invert();
					_Q_.Compress();
					// Qt.write(std::cout<< "Qtt = ") << std::endl;
				}

				if (_P_.isIdentity() && _Q_.isIdentity()) { /*  B and C are submatrices */
					_C_ = ( BlasMatrix<Integer>(A,0,0,_r_,_r_));
					_B_ = (BlasMatrix<Integer>(A,0,_r_,_r_,_s_));
				}
				else { /*  B and C should be extracted */
					_C_ = ( BlasMatrix<Integer>(_r_,_r_)) ;
					_B_ = ( BlasMatrix<Integer>(_r_,_s_)) ;
				}
				if (!_P_.isIdentity() || !_Q_.isIdentity()) {
					StackMat<BlasMatrix<Integer> >  QAP(_C_,_B_) ; // A = (C|B\?|?)
					extractMatrix(QAP,A,_Q_,_P_) ;
				}


			}

		}
		void getDims(size_t & r, size_t & s)
		{
			r = _r_ ;
			s = _s_ ;
		}
		void getPermuts(BlasPermutation<size_t> &P,BlasPermutation<size_t> &Q)
		{
			P = _P_ ;
			Q = _Q_ ;
		}
	} ;

	template<class Field>
	class NullspaceDomain {
		typedef typename Field::Element Element;
		BlasMatrix<Integer>    & _A_ ;
		// double                   _p_ ;
		// unsigned int             _n_ ;
		/*  */
		// BlasMatrix<Integer> & _Ared_ ;
		// NullspaceRep<Field>    & _N_ ;
		size_t                    _r_ ;
		size_t                    _s_ ;
		ReduceForm                _R_ ;

	public :
		ModularNullspace(BlasMatrix<Integer> & A) :
			_A_(A)
		{ }
		ModularNullspace(const BlasMatrix<Integer> & A) :
			_A_(A)
		{ }
		void initialize( )
		{
			_R_(A);
			// _n_ = n ;
			// _p_ = p ;
		}
		void reduce()
		{
			// Integer pn = pow(Integer(_p_),long(_n_));
			// Modular<Integer> Zpn(pn);
			// BlasMatrix<Integer> Ared(Zpn,_A_);
			// _Ared_ = Ared ;
		}

		void compute()
		{
			// RightDixon(_Ared_,_p_,_n_,_N_);
		}

#if 0
		NullspaceRep<Field> & getNullspace(NullspaceRep<Field> & Z)
		{
			return Z = _N_ ;
		}
#endif
	} ;
#endif

	/*! Space efficient representation of a nullspace in Hermite form (?).
	 * It is represented as \c P_.StackMatrix(Z_,m_*Id).
	 */
	template<class Field>
	class NullspaceRep {
		typedef typename Field::Element Element ;
	protected :
		DenseMatrix<Element>     & Z_ ;
		BlasPermutation<size_t>  & P_ ;
		Element                  & m_ ;
	public :
		NullspaceRep(DenseMatrix<Element> & Z,
			     BlasPermutation<size_t> & P,
			     Element & mid) :
			Z_(Z), P_(P), m_(mid)
		{
			if (!P.isIdentity()){
				linbox_check(P.getPointer() == P_.getPointer());
			}
		}

		NullspaceRep(DenseMatrix<Element> & Z,
			     BlasPermutation<size_t> & P) :
			Z_(Z), P_(P)
		{
			m_ = Element(-1);
		}

		NullspaceRep(NullspaceRep<Field> & N) :
			Z_(N.Z_),P_(N.P_),m_(N.m_)
		{}

		DenseMatrix<Element> getNullspace()
		{
			BlasMatrix<Element> N(Z_.coldim()+Z_.rowdim(),Z_.coldim());
			for (size_t i = 0 ; i < Z_.rowdim(); ++i)
				for (size_t j = 0 ; j < Z_.coldim(); ++j)
					N.setEntry(i,j,Z_.getEntry(i,j));

			for (size_t i = 0; i < Z_.coldim() ; ++i)
				N.setEntry(i+Z_.rowdim(),i,m_);
			Field F;
			BlasMatrixDomain<Field> BMD(F) ;
			BMD.mulin_right (P_,N);
		}

		NullspaceRep<Field> & operator=(NullspaceRep<Field> & N)
		{
			Z_(N.Z_);
			P_(N.P_);
			m_(N.m_);
			return *this;
		}

		/*!
		 * @pre Z_ has the form StackMatrix(Z,0*Id).
		 * @warning \c refNullspace() breaks the \c NullspaceRep
		 * compact representation and \p Z_ gets the full nullspace.
		 */
		DenseMatrix<Element> & refNullspace()
		{
			size_t r = Z_.rowdim()-Z_.coldim() ;
			for (size_t i = 0; i < Z_.coldim() ; ++i)
				Z_.setEntry(i+r,i,m_);
			Field F;
			BlasMatrixDomain<Field> BMD(F) ;
			BMD.mulin_right (P_,Z_);
		}

	} ;

	void PrepareNullspace(const BlasMatrix<Integer> & A,
			      BlasMatrix<Integer> & C,
			      BlasMatrix<Integer> & B,
			      BlasPermutation<size_t> & P,
			      BlasPermutation<size_t> & Qt)
	{

	}

	// write a p-adic matrix to o
	template<class T>
	std::ostream & write(std::ostream & o,
			     const std::vector<BlasMatrix<T > > & M, const size_t p,
			     const size_t subsize,
			     bool maple=false, const char * name = "A", const char * tmp_name = "ip")

	{
		Integer P = 1 ;
		if (!subsize)
			return o ;
		if (subsize==1) {
			if (!maple)
				M[0].write(o << name << ": ",maple);
			else {
#if 0
				o << "Mid := Matrix(" << M[0].rowdim() << ',' << M[0].coldim() << ",(i,j)->if (i=j+" << M[0].rowdim()-M[0].coldim() << ") then -1 else 0 end if ): " << std::endl ;
#endif
				o << name << ":= evalm(" ;
				M[0].write(o,maple)  ;
#if 0
				std::cout << "+ Mid)";
#endif

			}
			return o ;
		}
		if (!maple) {
			for (size_t i = 0 ; i < subsize ; ++i){
				o << "order " << P << ":" << std::endl;
				M[i].write(o,maple) ;
				o << std::endl;
				P *= p ;
			}
			return o ;
		}
		else {
			for (size_t i = 0 ; i < subsize ; ++i){
				M[i].write(o<<name<<tmp_name<<"["<<i+1<<"] := ",true)<<" : ";
			}
			o << std::endl;
#if 0
			o << "Mid := Matrix(" << M[0].rowdim() << ',' << M[0].coldim() << ",(i,j)->if (i=j+" << M[0].rowdim()-M[0].coldim() << ") then -1 else 0 end if ): " << std::endl ;
			o << name << ":= evalm(add((" <<name << tmp_name << "[i]+Mid)*"<<p<<"^(i-1),i=1.." << subsize << "))" << std::endl;
#endif
			o << name << ":= evalm(add((" <<name << tmp_name << "[i])*"<<p<<"^(i-1),i=1.." << subsize << "))";
			return o ;

		}
	}



	// P = P_0 + P_1 p +...+ P_k p^k
	template<class DMatrix, class IMatrix>
	std::vector<DMatrix> & Padize(std::vector<DMatrix> & P
				      , const IMatrix & M
				      , const double p)
	{
		size_t m = M.rowdim();
		size_t n = M.coldim();
		// unsigned long sz     = 0 ;
		// unsigned long tmp_sz = 0 ;
		/*  get max size */
		unsigned long sz ;
		MaxElementSize(sz,M);
#if 0 /*  faster ? */
		Integer big ;
		MaxElementAbs(big,M);
#endif
		size_t big = ceil((double) sz/(double) logtwo(p));
		// std::cout << "#largest element of size " << sz << " divided into " << big << " chunks" << std::endl;
		P.resize(big);
		const IMatrix Zero(m,n);
		for (size_t i = 0 ; i < big ; ++i)
			P[i] = Zero ;

		Integer dvd ;
		Integer quo ;

		typedef Modular<double> Field;
		typedef Field::Element Element;
		long rem ;

		for (size_t i = 0 ; i < M.rowdim() ; ++i)
			for (size_t j = 0 ; j < M.coldim() ; ++j) {
				size_t k   = 0;
				dvd     = M.getEntry(i,j);
				while(1) {
					assert(k<big);
					Integer::divmod(quo,rem,dvd,p);
#ifdef DEBUG
					if (dvd != Integer(p)*quo+Integer(rem))
						std::cout << dvd << '=' << p << " * (" << quo << ") + (" << rem << ");" << std::endl;
#endif
					linbox_check(dvd == Integer(p)*quo+Integer(rem));
					linbox_check(rem>=0 && (Element)rem<p);
					// linbox_check(quo>=0);
					P[k].setEntry(i,j,(Element) rem) ;
					++k;
					if (abs(quo) < (Integer)p) {
						assert(k<big);
						P[k].setEntry(i,j,(Element)quo) ;
						break;
					}
					dvd = quo ;
				}
			}

#if 0 /* Alternative */
		IMatrix G(M);
		for (size_t i = 0 ; i < big ;++i) {
			for (Iter Git = G.rawBegin(); Git != G.rawEnd();  ++Git){
				if (*Git){
					Integer::divmod(P[i],*Git,*Git,p);
				}
			}
		}
#endif

		// write(std::cout, P,p,true,"P_decomp") <<  std::endl ;
		return P ;


	}

	template<class DMatrix, class IMatrix>
	std::vector<DMatrix> & Padize(std::vector<DMatrix> & P
				      , const IMatrix & M
				      , const double p
				      , const size_t d)
	{
		size_t m = M.rowdim();
		size_t n = M.coldim();
		linbox_check(P.size() == d);
		linbox_check(P[0].rowdim()==m);
		linbox_check(P[0].coldim()==n);

		Integer dvd ;
		Integer quo ;

		typedef Modular<double> Field;
		typedef Field::Element Element;
		long rem ;

		for (size_t i = 0 ; i < m ; ++i)
			for (size_t j = 0 ; j < n ; ++j) {
				size_t k   = 0;
				dvd     = M.getEntry(i,j);
				while(1) {
					Integer::divmod(quo,rem,dvd,p);
#ifdef DEBUG
					if (dvd != Integer(p)*quo+Integer(rem))
						std::cout << dvd << '=' << p << " * (" << quo << ") + (" << rem << ");" << std::endl;
#endif
					linbox_check(dvd == Integer(p)*quo+Integer(rem));
					linbox_check(rem>=0 && (Element)rem<p);
					// linbox_check(quo>=0);
					P[k].setEntry(i,j,(Element) rem) ;
					++k;
					if (abs(quo) < (Integer)p) {
						if(k<d)
							P[k].setEntry(i,j,(Element)quo) ;
						break;
					}
					if (k==d)
						break;
					dvd = quo ;
				}
			}

#if 0 /* Alternative */
		IMatrix G(M);
		for (size_t i = 0 ; i < big ;++i) {
			for (Iter Git = G.rawBegin(); Git != G.rawEnd();  ++Git){
				if (*Git){
					Integer::divmod(P[i],*Git,*Git,p);
				}
			}
		}
#endif

		// write(std::cout, P,p,true,"P_decomp") <<  std::endl ;
		return P ;


	}

	template<class DMatrix, class IMatrix>
	IMatrix & Dadize( IMatrix & M,
			  const std::vector<DMatrix> & P,
			  const double p  )
	{
		// typedef typename IMatrix::RawIterator Iter ;
		/*  get max size */
		linbox_check(P.size());
		size_t lig = P[0].rowdim();
		size_t col = P[0].coldim();
		linbox_check(lig<=M.rowdim());
		linbox_check(col<=M.coldim());
		for (size_t i = 0 ; i < lig ; ++i)
			for (size_t j = 0 ; j < col ; ++j)
				Integer::addin( M.refEntry(i,j) , Integer(P[0].getEntry(i,j)) ) ;

		Integer p_k = p ;

		for (size_t k = 1 ; k < P.size() ; ++k) {
			for (size_t i = 0 ; i < lig ; ++i){
				for (size_t j = 0 ; j < col ; ++j)
					Integer::axpyin(M.refEntry(i,j) , Integer(P[k].getEntry(i,j)) , p_k) ;

			}
			p_k *= p ;
		}
		// M.write(std::cout << "M_rebuild :=" ,true) << ';' << std::endl;

		return M ;


	}

}

namespace LinBox
{ /* IML nullspace */

	/*  IntegerNS */
	/*!@internal NO DOC
	 *
	 * @param n
	 * @param m
	 * @param A
	 * @param mp_N_pass
	 * @todo to be faster, consider the transpose of NN (vectors will be rows)
	 */
	size_t
	right_IntegerNS(
			BlasMatrix<Integer>         & N
			// , std::vector<Integer>      & d
			, const BlasMatrix<Integer> & A
			// , bool                        check = 0
		       )
	{
		size_t lA ;
		MaxElementSize(lA,A) ;
		std::vector<Integer> d ;

		size_t m = A.rowdim();
		size_t n = A.coldim();
		size_t r = 0  ; // rank
		size_t s ; // ker dim, the result.

		//typedef double Element ;
		typedef Modular<double> Field;     // dans lequel on réduit A
		typedef Field::Element  Element ;

		PID_integer ZZ ;
		MatrixDomain<PID_integer>        ZMD(ZZ);
		BlasMatrixDomain< PID_integer > BZMD(ZZ) ;

		size_t essai = _LB_MAX_NS_TRY ;

		while (1) { /* main loop */
			BlasPermutation<size_t> P (n);
			BlasPermutation<size_t> Qt(m);
			// ReduceForm RR(A);
			// RR.getDims(r,s);
			// RR.getPermuts(P,Qt);
			// RR.getCB(C,B);

			//!@todo use MatrixHom::map ?
			{ /* LU mod p  */
				RandomPrimeIter PrimeGen(19) ;
				Element p = (Element) PrimeGen.random_between(15);
				Field F(p);
				BlasMatrix<Element> Ap(m,n);
				MatrixHom::map(Ap,A,F);
				// Ap.write(std::cout << "A mod " << p << " :=",true) << std::endl;
				LQUPMatrix<Field> LUp(F,Ap,P,Qt);
				r = LUp.getrank();

			}

			s = n - r; // taille du noyal.
			d.resize(s);


			if (s == 0) { // A était inversible !
				N= BlasMatrix<Integer>(0,0);
				break;
			}
			else if (r == 0) { // A est nulle ?
				bool nulle = ZMD.isZero(A) ;
				if (!nulle) {
					--essai;
					std::cout << "Bad prime (very bad)" << std::endl;
					if (!essai)
						throw(LinBoxFailure(__func__,__FILE__,__LINE__,
								    "could not find a nullspace..."));
					continue; // A non nulle, pas de bol hein ?
				}
				N = BlasMatrix<Integer>(n,n);
				Integer one(1);
				for (size_t i = 0 ; i < n ; ++i){
					N.setEntry(i,i,one);
				}
				break ;
			}
			else { // A a un noyau non trivial !

				BlasMatrix<Integer> * C = NULL ;
				BlasMatrix<Integer> * B = NULL ;

				P.Compress();
				Qt.Compress();

				if (!Qt.isIdentity()) {
					Qt.Invert();
					Qt.Compress();
				}

				if (P.isIdentity() && Qt.isIdentity()) {
					C = (new BlasMatrix<Integer>(A,0,0,r,r));
					B = (new BlasMatrix<Integer>(A,0,r,r,s));
				}
				else {
					C = (new BlasMatrix<Integer>(r,r)) ;
					B = (new BlasMatrix<Integer>(r,s)) ;
				}

				BlasMatrix<Integer> NN(n,s) ; // le noyal !!

				StackMat<BlasMatrix<Integer> >  QAP(*C,*B) ; // A = (C|B\?|?)
				if (!P.isIdentity() || !Qt.isIdentity()) {
					extractMatrix(QAP,A,Qt,P) ;
				}


				Timer Solver ; Solver.clear();
				Solver.start();

#ifdef __LINBOX_USE_IML
				mpz_t * C_iml = reinterpret_cast<mpz_t*>(C->getWritePointer());
				mpz_t * B_iml,*N_iml, D_iml;
				mpz_init(D_iml);
#ifndef __LINBOX_CUTS_HAIR_IN_FOUR
				B_iml = reinterpret_cast<mpz_t*>(B->getWritePointer());
				N_iml = reinterpret_cast<mpz_t*>(NN.getWritePointer()) ;
				nonsingSolvLlhsMM(RightSolu, r, s, C_iml, B_iml, N_iml, D_iml);
				Solver.stop(); std::cout << "total solver : " << Solver << std::endl;
				Integer di ;
				mpz_set(di.get_mpz(),D_iml);
				Integer::negin(di) ;

				for (size_t i = 0 ; i < s ; ++i) {
					NN.setEntry( r+i , i , di) ;
				}
				mpz_clear(D_iml);
#else
				/*  B^t */
				BlasMatrix<Integer> Bt(s,r) ;
				for (size_t i = 0 ; i < s ; ++i)
					for (size_t j = 0 ; j < r ; ++j)
						Bt.setEntry( i,j,B->getEntry(j,i) );
				BlasMatrix<Integer> Nt(s,r) ;
				B_iml = reinterpret_cast<mpz_t*>(Bt.getWritePointer());
				N_iml = reinterpret_cast<mpz_t*>(Nt.getWritePointer()) ;
				for (size_t j = 0 ; j < s ; ++j) {
					nonsingSolvLlhsMM(RightSolu, r, 1, C_iml, B_iml, N_iml, D_iml);
					B_iml += r;
					N_iml += r;
					mpz_set(d[j].get_mpz(),D_iml);
				}
				for (size_t i = 0 ; i < r ; ++i)
					for (size_t j = 0 ; j < s ; ++j)
						NN.setEntry( i,j,Nt.getEntry(j,i) );

				for (size_t i = 0 ; i < s ; ++i) {
					Integer::negin(d[i]) ;
					NN.setEntry( r+i , i , d[i]) ;
				}

				Solver.stop(); std::cout << "total solver : " << Solver << std::endl;
#endif
#else
				// d.CxN=B
				RationalSolver<PID_integer,Modular<double>,RandomPrimeIterator,DixonTraits> QD ;
				RationalSolver<PID_integer,Modular<double>,RandomPrimeIterator,NumericalTraits> QS ;
				if (lA<50)
					QS.hasSmallCoeffs(); //!@bug B two...

				BlasBlackbox<PID_integer> CC (ZZ,QAP.first);
#if 0
				Timer Solver_in ; Solver_in.clear();
				Timer Chrono ;
#endif

				bool recontructed = true ;
				for (size_t j = 0 ; j < s ; ++j) {
					std::vector<Integer> v_n(r);
					std::vector<Integer> v_b(r);
					for (size_t i = 0 ; i < r ; ++i)
						v_b[i] = B->getEntry(i,j) ;
					// QQ.solveNonsingular(v_n,d,QAP.first,v_b) ; //! @todo
					enum SolverReturnStatus SRS ;
#if 0
					Chrono.clear(); Chrono.start();
#endif
					if (lA<50) {
						// std::cout << "trying fast solve" << std::endl;
						SRS = QS.solve(v_n,d[j],CC,v_b) ; //! @todo
						if (SRS != SS_OK){
							// std::cout << "failed, falling back to standard" << std::endl;
							SRS = QD.solveNonsingular(v_n,d[j],CC,v_b,false,5) ; //! @todo
						}
					}
					else {
						SRS = QD.solveNonsingular(v_n,d[j],CC,v_b,false,5) ; //! @todo
					}
#if 0
					Chrono.stop(); Solver_in += Chrono ;
#endif
#if 0 /*  debug */
					std::vector<Integer> v_r(r);
					CC.apply(v_r,v_n);
					for (size_t i = 0 ; i < r ; ++i)
						if (v_b[i]*d[j] != v_r[i])
							std::cout << 'X' ;
						else
							std::cout << 'O' ;
					std::cout << '\n' ;
#endif
					if (SRS != SS_OK){
						--essai;
						std::cout << "Recoonstruction failed" << std::endl;
						if (!essai)
							throw(LinBoxFailure(__func__,__FILE__,__LINE__,"could not find a nullspace..."));

						recontructed = false ;
						break;

					}
					for (size_t i = 0 ; i < r ; ++i)
						NN.setEntry(i,j,v_n[i])  ;

				}
#if 0
				std::cout << "only solver : " << Solver_in << std::endl;
				std::cout << "per iter : " << Solver_in.realtime()/s << std::endl;
				Solver.stop(); std::cout << "total solver : " << Solver << std::endl;
#endif
				if (!recontructed){
					continue ;
				}
				delete B ;
				delete C ;


				// on rajoute l'identité
				for (size_t i = 0 ; i < s ; ++i) {
					Integer::negin(d[i]) ;
					NN.setEntry( r+i , i , d[i]) ;
				}
#endif


				/* on applique P^{-1} */
				TransposedBlasMatrix<BlasPermutation<size_t > > PP(P);
				BZMD.mulin_right (PP,NN);


				{/*  Checking with other rows */
					size_t reste = m-r;
					const size_t t = std::log(reste)+1;
					BlasMatrix<Integer>Ax(t,n);

					for (int i=t ; i-- ;) {
						//!@todo comment faire sans copie ? créer deux vecteurs 'statiques' et jouer aux pointeurs ?
						int i_q = Qt[r+Integer::random_lessthan(Integer(reste))] ;
						linbox_check(!(i_q<0) && i_q<(int)m);
						for ( size_t j = 0 ; j < A.coldim() ; ++j)
							Ax.setEntry(i,j, A.getEntry(i_q,j));
					}

					RandomPrimeIter PrimeGen(19) ;
					Element p = (Element) PrimeGen.random_between(15);
					Field F(p);
					BlasMatrix<Element>     Axp(Ax,F);
					BlasMatrix<Element>      Np(NN,F);
					BlasMatrix<Element>       Rp(t,s);
					BlasMatrixDomain<Field>   BFMD(F);
					MatrixDomain<Field>        FMD(F);

					BFMD.mul(Rp,Axp,Np);
					if(!FMD.isZero(Rp)){
#if 0 /*  debug */
						std::cout << "p := " << p << ";" << std::endl;
						Axp.write(std::cout << "Axp :=  ",true) << ';' << std::endl;
						Np.write(std::cout << "Np :=  ",true) << ';' << std::endl;
						Rp.write(std::cout << "Rp :=  ",true) << ';' << std::endl;
#endif
						--essai ;
						std::cout << "trying again..." << std::endl;
						if (!essai)
							throw(LinBoxFailure(__func__,__FILE__,__LINE__,"could not find a nullspace..."));

						continue;
					}
				}

				N = NN ;
				break;
			}
		}
		return s;

	}

	size_t IntegerNS(
			 BlasMatrix<Integer>         & N
			 // , std::vector<Integer>      & d
			 , const BlasMatrix<Integer> & A
			 , const LinBoxTag::Side Side
			)
	{
		if (Side == LinBoxTag::Right)
			return right_IntegerNS(N/*,d*/,A);
		else
			return -1 ;
	}

}

namespace LinBox
{ /*  mod p^n */

	int Dixon(std::vector<BlasMatrix<double> > & N,
		  const BlasMatrix<Integer> & Cz, const BlasMatrix<Integer> &Bz,
		  double p, unsigned int d)
	{
		typedef Modular<double> Field;     // dans lequel on réduit A
		typedef Field::Element  Element ;

		Field F(p);
		PID_integer ZZ ;

		BlasMatrixDomain<Field> PMD(F);
		MatrixDomain<PID_integer> ZMD(ZZ);


		size_t r = Cz.rowdim();
		size_t s = Bz.coldim();

		Integer pd = pow(Integer(p),long(d));
		Modular<Integer> Zpd(pd);
		BlasMatrix<Integer> C(Cz,Zpd);
		BlasMatrix<Integer> B(Bz,Zpd);



		BlasMatrix<Element> Cp(C,F); // Cp = C mod p
		BlasMatrix<Element> Bp(B,F); // Bp = B mod p
		BlasMatrix<Element> Cinv(r,r);
		int null = 0 ;
		PMD.inv(Cinv,Cp,null);
		if (null){
			std::cout << "failed to invert C" << std::endl;
			return -2 ;
		}

		BlasMatrix<Integer> R(r,s);
		Integer ZZone(1);
		Integer ZZmin(-1);

		{ /* main loop  */
			BlasMatrix<Integer> AN(r,s);
			for (size_t i = 0 ; i < (size_t)d ; ++i){
				BlasMatrix<Element> Ni(N[i]);
				/*  A N_i = R mod p */
				MatrixHom::map(Ni,R,F);
				PMD.addin(Ni,Bp) ;
				PMD.mulin_right(Cinv,Ni);
				/*  update R */
				if (i+1<(size_t)d) {
					//!@todo ZMD.muladdin !!!!!! qui remplace tout ça :
					ZMD.mul(AN,C, Ni);
					ZMD.subin(R,AN);
					ZMD.addin(R,B);
					BlasMatrix<Integer>::RawIterator it = R.rawBegin();
					for ( ; it != R.rawEnd() ; ++it ) {
						*it /= p ;
					}
				}
				linbox_check(Ni.getEntry(0,0) == N[i].getEntry(0,0));

			}
		}

		return 0 ;

	}

	int MoenckCarter(std::vector<BlasMatrix<double> > & N,
			 const BlasMatrix<Integer> & Cz, const BlasMatrix<Integer> &Bz,
			 double p, unsigned int d)
	{
		typedef Modular<double> Field;     // dans lequel on réduit A
		typedef Field::Element  Element ;

		Field F(p);
		PID_integer ZZ ;

		BlasMatrixDomain<Field> PMD(F);
		MatrixDomain<PID_integer> ZMD(ZZ);

		Integer pd = pow(Integer(p),long(d));
		Modular<Integer> Zpd(pd);
		BlasMatrix<Integer> C(Cz,Zpd);
		BlasMatrix<Integer> B(Bz,Zpd);

		size_t r = C.rowdim();
		size_t s = B.coldim();

		BlasMatrix<Element> Cp(C,F); // Cp = C mod p
		BlasMatrix<Element> Bp(B,F); // Bp = B mod p
		BlasMatrix<Element> Cinv(r,r);
		int null = 0 ;
		PMD.inv(Cinv,Cp,null);

		if (null){
			std::cout << "failed to invert C" << std::endl;
			return -2 ;
		}


		// BlasMatrix<Integer> R(r,s); // Delta
		Integer ZZone(1);
		Integer ZZmin(-1);


		{ /* main loop  */
			BlasMatrix<Integer> AX(r,s);
			// BlasMatrix<Element> AX_temp(r,s);
			BlasMatrix<Integer> AY_temp(r,s);
			BlasMatrix<Element> X(r,s);
			BlasMatrix<Integer> Bi(r,s);

			Integer pi = p ;
			Integer un = 1 ;
			Integer pim1 = 1 ;

			for (size_t k = 0 ; k < (size_t)d ; ++k){
				BlasMatrix<Element> Nk(N[k]);
				if (k) {
					BlasMatrix<Element> Nkm1(N[k-1]);
					/* AX +=  p^{i-1}*A X */
					ZMD.mul(AY_temp,C,Nkm1) ; // Cp2 . X
					ZMD.subin(AY_temp,B);// - Bp2

					for (size_t i = 0 ; i < r ; ++i)
						for (size_t j = 0 ; j < s ; ++j) {
							// AX.refEntry(i,j) -= Bp2.getEntry(i,j);
							Integer::axpyin(AX.refEntry(i,j),pim1,
									AY_temp.getEntry(i,j));  // + AX
						}
					pim1 *= p ;
					// std::cout << "# " ;
					for (size_t i = 0 ; i < r ; ++i)
						for (size_t j = 0 ; j < s ; ++j) {
							// if (AX.getEntry(i,j) % pim1)
							// std::cout << "f "  ;
							// else
							// std::cout << "v "  ;
							F.init(Nk.refEntry(i,j),-AX.getEntry(i,j) / pim1);
						}
				}
				// std::cout << std::endl;
				/* Y_i = C W mod p^2 */
				PMD.addin(Nk,Bp);
				PMD.mulin_right(Cinv,Nk);

				linbox_check(Nk.getEntry(0,0) == N[k].getEntry(0,0));

			}
		}
		// write(std::cout,NN,p,NN.size(),true,"Z") << ';' << std::endl;

		return 0 ;

	}

	//!@todo make this class "follow" \c l.
	template<class Field>
	class Moduleur {
		typedef typename Field::Element Element ;
		typedef UnparametricField<double> DoubleDomain ;
	private :
		std::vector<BlasMatrix<Element > > & _M_ ;
		std::vector<BlasMatrix<Element > > & _C_ ;
		size_t _d_ ;
		size_t _r_ ;
		size_t _s_ ;
		long   _p_ ;
		double _more_food ;
		DoubleDomain _D_ ;
		BlasMatrixDomain<DoubleDomain> _BMD_;


		void digest(int l)
		{
			// std::cout << "# cleaning at " << l << std::endl;
			// write(std::cout,_M_,_p_,_d_,true,"before") << ";" << std::endl;
			ldiv_t qr ;
			for(size_t i = 0 ; i < _r_ ; ++i)
				for(size_t j = 0 ; j < _s_ ; ++j) {
					size_t k = l ;
					while ( k < _d_ && !(_M_[k].getEntry(i,j)<_p_) ) {
						qr = std::div((long)_M_[k].getEntry(i,j),(long)_p_);
						_M_[k].setEntry(i,j,qr.rem);
						if (k+1<_d_)
							_M_[k+1].refEntry(i,j) += qr.quot ;
						++k;
					}
					// if (k>(size_t)l)
					// std::cout << k-l << " " ;

				}
			// std::cout << std::endl;
			// write(std::cout,_M_,_p_,_d_,true,"after") << ";" << std::endl;
			return;
		}

	public :
		Moduleur(const Field & F
			 , std::vector<BlasMatrix<Element> > & Bdic
			 , std::vector<BlasMatrix<Element> > & Cdic
			) :
			_M_(Bdic),_C_(Cdic),_BMD_(_D_)
		{
			_d_ = Bdic.size();
			linbox_check(_d_);
			_r_ = Bdic[0].rowdim();
			_s_ = Bdic[0].coldim();
			_p_ = F.characteristic();
			_more_food = F.AccBound( 1 );
			linbox_check(_more_food > _r_ );
			_more_food = int(_more_food)/ _r_ ;

			// DoubleDomain D;
			// _BMD_(D);
			// write(std::cout,_M_,_p_,_d_,true,"mBadic") << ";" << std::endl;
			// write(std::cout,_C_,_p_,_d_,true," Cadic") << ";" << std::endl;

		}

		void feed( const std::vector<BlasMatrix<Element> > & N,
			   int l)
		{
			linbox_check(l);

			// std::cout << "# before feeding " << l << std::endl;
			// write(std::cout,_M_,_p_,min((size_t)l+1,_d_),true,"Videu") << ";" << std::endl;

			linbox_check((size_t)l < _d_);
			/*  on ajoute la suite */
			size_t food = _more_food ;
			// BlasMatrix<Element> *Acc0 = &(_M_[0]) ;
			_BMD_.axpyin(_M_[l-1],_C_[0],N[l-1]);
			digest(l-1);
			for (size_t j = 0 ; j < (size_t)l ; ++j) {
				// std::cout << N[j].getEntry(0,0) << " = " ;
				_BMD_.axpyin(_M_[l],_C_[l-j],N[j]);
				// std::cout << N[j].getEntry(0,0) << std::endl ;
				--food;
				if (!food){
					digest(l);
					food = _more_food ;
				}
			}
			digest(l);
			// std::cout << "# after feeding " << l << std::endl;
			// write(std::cout,_M_,_p_,min((size_t)l+1,_d_),true,"Plein") << ";" << std::endl;
			return;
		}

		void pee(BlasMatrix<Element > & Delta, int l)
		{
			Delta =(BlasMatrix<Element > & ) _M_[l] ;
		}
	};

	int Noix(std::vector<BlasMatrix<double> > & N,
		 const BlasMatrix<Integer> & C, const BlasMatrix<Integer> &B,
		 double p, unsigned int d)
	{
		// std::cout << "p,d:=" << p << ',' << d << ':' <<std::endl;
		typedef Modular<double> Field;     // dans lequel on réduit A
		typedef Field::Element  Element ;

		Field F(p);
		PID_integer ZZ ;

		BlasMatrixDomain<Field> PMD(F);
		MatrixDomain<PID_integer> ZMD(ZZ);

		size_t r = C.rowdim();
		size_t s = B.coldim();


		// B.write(std::cout << "B :=", true) << ';' << std::endl;
		// C.write(std::cout << "C :=", true) << ';' << std::endl;

		/*  images of C and B in Z/(p^d)Z */
		Integer pd = pow(Integer(p),long(d));
		Modular<Integer> Zpd(pd);
		const BlasMatrix<Element> Ze(r,s);
		const BlasMatrix<Element> Ro(r,r);
		std::vector<BlasMatrix<Element> > Badic(d,Ze);
		std::vector<BlasMatrix<Element> > Cadic(d,Ro);
		{
			/*  Bdic *= (p^d-1)/(p-1) */
			long p1 = p-1 ;
			Integer q,r ;
			//ppd1 = (p^d-1)/(p-1)
			Integer ppd1 = pow(Integer(p),long(d));
			Integer u ;
			inv(u,p1,ppd1);

			BlasMatrix<Integer> Cpd(C,Zpd);
			BlasMatrix<Integer> Bpd(B.rowdim(),B.coldim());
			for (size_t i = 0 ; i< Bpd.rowdim(); ++i) // -B
				for (size_t j = 0 ; j< Bpd.coldim(); ++j) {
					Zpd.init(Bpd.refEntry(i,j),B.getEntry(i,j));
					Integer::mulin(Bpd.refEntry(i,j),u);
					// Integer::divmod(q,r,Bpd.getEntry(i,j),p1);
					// #ifdef DEBUG
					// if ( Bpd.getEntry(i,j)!= Integer(p1)*q+r)
					// std::cout << "!!!!" << Bpd.getEntry(i,j) << '=' << p1 << " * (" << q << ") + (" << r << ");" << std::endl;
					// #endif

					// Zpd.init(Bpd.refEntry(i,j),Integer(r*ppd1-q));
				}
			Padize(Badic,Bpd,p,d);
			Padize(Cadic,Cpd,p,d);
		}
		// write(std::cout << "B-adic := ",Badic,p,d,true,"Badic") << ";" << std::endl;
		// write(std::cout << "C-adic := ",Cadic,p,d,true,"Cadic") << ";" << std::endl;

		BlasMatrix<Element> Cp(C,F); // Cp = C mod p
		BlasMatrix<Element> Bp(B,F); // Bp = B mod p
		// Bp.write(std::cout << "Bp :=",true) << ';' << std::endl;

		/*  inverse C mod p */
		BlasMatrix<Element> Cinv(r,r);
		int null = 0 ;
		PMD.inv(Cinv,Cp,null);
		if (null){
			std::cout << "failed to invert C" << std::endl;
			return -2 ;
		}

		Integer ZZone(1);
		Integer ZZmin(-1);

		Moduleur<Field> Delta(F,Badic,Cadic);

		{ /* main loop  */
			BlasMatrix<Element> R(r,s);
			BlasMatrix<Integer> AX(r,s);
			// BlasMatrix<Element> AX_temp(r,s);
			BlasMatrix<Integer> AY_temp(r,s);
			BlasMatrix<Element> X(r,s);
			BlasMatrix<Integer> Bi(r,s);
			MatrixDomain<Field> MD(F);

			Integer pi = p ;
			Integer un = 1 ;
			Integer pim1 = 1 ;

			for (size_t k = 0 ; k < (size_t)d ; ++k){
				// std::cout << "# *** stage " << k << " started ***" << std::endl;
#if 1
				BlasMatrix<Element> Nk((const BlasMatrix<Element>&)N[k]);
#endif
				if (k) {
					Delta.feed(N,k);
					Delta.pee(N[k],k);
					// linbox_check(Nk.getEntry(0,0) == N[k].getEntry(0,0));
					MD.negin(N[k]);
					// N[k].write(std::cout << "Delta_" << k << " :=",true) << ';' << std::endl;
				}
#if 1
				if (k) {
					BlasMatrix<Element> Nkm1(N[k-1]);
					/* AX +=  p^{i-1}*A X */
					ZMD.mul(AY_temp,C,Nkm1) ; // Cp2 . X
					ZMD.subin(AY_temp,B);// - Bp2

					for (size_t i = 0 ; i < r ; ++i)
						for (size_t j = 0 ; j < s ; ++j) {
							// AX.refEntry(i,j) -= Bp2.getEntry(i,j);
							Integer::axpyin(AX.refEntry(i,j),pim1,
									AY_temp.getEntry(i,j));  // + AX
						}
					pim1 *= p ;
					// AX.write(std::cout << "Reference_" << k << ":=",true) << ';' << std::endl;
					// std::cout << "# " ;
					for (size_t i = 0 ; i < r ; ++i)
						for (size_t j = 0 ; j < s ; ++j) {
							// if (AX.getEntry(i,j) % pim1)
							// std::cout << "f "  ;
							// else
							// std::cout << "v "  ;
							F.init(R.refEntry(i,j),-AX.getEntry(i,j) / pim1);
						}
				}
				// std::cout << std::endl;
				// R.write(std::cout << "Refer_" << k << ":=", true) << ';' << std::endl;
#endif
				if (!k)
					PMD.addin(N[k],Bp);
				PMD.mulin_right(Cinv,N[k]);
				// N[k].write(std::cout << "Z_" << k << " :=",true) << ';' << std::endl;
				// write(std::cout ,N,p,k+1,true,"ZZ") << ";" << std::endl;
				// std::cout << "# *** stage " << k << " ended ***" << std::endl;
				// linbox_check(Nk.getEntry(0,0) == N[k].getEntry(0,0));

			}
		}

		return 0 ;

	}


	enum NullSpaceMethod {
		Dixon_,
		MoenckCarter_,
		Noix_
	} ;

	template<class IMatrix, class DMatrix>
	int RightPadicNullspace(std::vector<DMatrix >& N, const IMatrix & C, const IMatrix & B,
			       	const double p, const int d,
				const enum NullSpaceMethod & meth)
	{
		int erreur = 0 ;
		switch(meth) {
		case (Dixon_) :
			{
				erreur = Dixon(N,C,B,p,d);
				break;
			}
		case (MoenckCarter_) :
			{
				erreur = MoenckCarter(N,C,B,p,d);
				break;
			}
		case (Noix_) :
			{
				erreur = Noix(N,C,B,p,d);
				break;
			}

		default :
			{
				erreur = -1 ;
				std::cout << "unknown P-adic method." << std::endl;
			}

		}
		return erreur ;
	}

	// A.N = 0 mod p^d
	template<class IMatrix>
	int RightPadicNullspace(IMatrix & N, const IMatrix & A,
				const double p, const int d,
			       	const NullSpaceMethod & meth)
	{
		size_t m = A.rowdim();
		size_t n = A.coldim();
		size_t r = 0  ;    // rank
		size_t s ;         // ker dim, the result.

		typedef Modular<double> Field;     // dans lequel on réduit A
		typedef Field::Element  Element ;

		PID_integer ZZ ;
		Field F(p);

		MatrixDomain<PID_integer> ZMD(ZZ);

		BlasPermutation<size_t> P (n);
		BlasPermutation<size_t> Qt(m);

		{ /* LU mod p  */
			//!@todo reuse this info !!
			BlasMatrix<Element> Ap(m,n);
			MatrixHom::map(Ap,A,F);
			LQUPMatrix<Field> LUp(F,Ap,P,Qt);
			r = LUp.getrank();
		}

		s = n - r; // taille du noyal.
		// std::cout << "ker dim is (probably) " << s << std::endl;
		// d.resize(s);

		if (s == 0) { // A était inversible !
			// std::cout << "invertible !" << std::endl;
			N= BlasMatrix<Integer>(0,0);
			/* A is Probably invertible. second chance ? */
			// std::cout << "not checking again" << std::endl;
			return s ;
		}
		else if (r == 0) { // A est nulle ?
			// std::cout << "nulle !" << std::endl;
			bool nulle = ZMD.isZero(A) ;
			if (!nulle) {
				std::cout << "error : bad p !" << std::endl;
				return -1 ; // fail
			}
			N = BlasMatrix<Integer>(n,n);
			Integer one(1);
			for (size_t i = 0 ; i < n ; ++i){
				N.setEntry(i,i,one);
			}
			// std::cout << "ok, kernel is Id" << std::endl;
			return s ;
		}
		else { // A a un noyau non trivial !
			// std::cout << "on a un noyal !" << std::endl;

			/* pre-allocating it */
			const BlasMatrix<Element > Zerooo(r,s); // add id ?
			std::vector<BlasMatrix<Element > > NN(d,Zerooo); // le noyal !!

			/*  extracting C/B */
			BlasMatrix<Integer> * C = NULL ;
			BlasMatrix<Integer> * B = NULL ;

			P.Compress();
			Qt.Compress();

			if (!Qt.isIdentity()) {
				Qt.Invert();
				Qt.Compress();
			}

			if (P.isIdentity() && Qt.isIdentity()) {
				C = (new BlasMatrix<Integer>(A,0,0,r,r));
				B = (new BlasMatrix<Integer>(A,0,r,r,s));
			}
			else {
				C = (new BlasMatrix<Integer>(r,r)) ;
				B = (new BlasMatrix<Integer>(r,s)) ;
			}

			StackMat<BlasMatrix<Integer> >  QAP(*C,*B) ; // A = (C|B\?|?)
			if (!P.isIdentity() || !Qt.isIdentity()) {
				extractMatrix(QAP,A,Qt,P) ;
			}

			/*  Nullspace mod p^d */
			int erreur =  RightPadicNullspace(NN,*C,*B,p,d,meth);
			if (erreur){
				std::cout << " *** P-adic solver has had errors *** " << std::endl;
			}

			/*  aggregating parts  */

			// write(std::cout << "Z := ",NN,p,true) << ";" << std::endl;
			BlasMatrix<Integer> K(n,s);
			N = Dadize(K,NN,p);
			Integer pn = Integer(1)-pow((Integer)p,(long)d);
			pn /= (Integer)(p-1);
			for (size_t i = 0; i < s; i++)
				N.setEntry( r+i , i , pn) ;
			BlasMatrixDomain<PID_integer> ZZDom(ZZ);
			TransposedBlasMatrix<BlasPermutation<size_t > > Pt(P);
			ZZDom.mulin_right (Pt,N);

			// N.write(std::cout << "N := ",true) << ';' << std::endl;


			/*-------------------*/
		}

		return s ;

	}

	template<class IMatrix>
	class RightNullspace {
		typedef Modular<double>  ModularField;
		typedef Modular<Integer> IntegerField;
		PID_integer ZZ ;
		IMatrix                     & _A_ ;
		IMatrix                     * _C_ ;
		IMatrix                     * _B_ ;
		BlasPermutation<size_t>     & _P_ ;
		BlasPermutation<size_t>    & _Qt_ ;
		NullspaceRep<IntegerField>  & _Z_ ;

	public:
		RightNullspace(const IMatrix & A) :
			_A_(A)
		{}

		void prepare()
		{
			// create _C_,...
		}
	};
}

namespace LinBox
{ /* p-adic + CRA */

#if 1 /*  bad LU mod p + CRA idea */
	template<class Field>
	struct NS_mod_pd {
		typedef typename Field::Element Element;
		// long int r_  ;
		// size_t   m_ ;
		// size_t   n_ ;
		BlasMatrix<Integer> & _C_ ;
		BlasMatrix<Integer> & _B_ ;
		size_t _d_ ;
		enum NullSpaceMethod _meth_ ;
		size_t _r_ ;
		size_t _s_ ;
		size_t _n_ ;

		NS_mod_pd (BlasMatrix<Integer> & C, BlasMatrix<Integer> & B, size_t d,
			   NullSpaceMethod  meth) :
			_C_(C), _B_(B), _d_(d), _meth_(meth),
			_r_(_C_.rowdim()), _s_(B.coldim()), _n_(_r_+_s_)
		{} ;

		BlasMatrix<Integer>& operator()(BlasMatrix<Integer>& N, const Field& F) const
		{
			Integer p = F.characteristic();
			std::cout << "P:=" << p << std::endl;
			double _p_ = naturallog(p)/_d_ ;
			_p_ = (long) exp(_p_);
			// std::cout << "prime :" << _p_ << std::endl;
			if (!(pow((integer)_p_,_d_)==p)) {
				_p_ += 1 ;
				// std::cout << "+1 :" << _p_ << std::endl;
			}
			if (!(pow((integer)_p_,_d_)==p)){
				_p_ -= 2 ;
				// std::cout << "-1 :" << _p_ << std::endl;
			}
			// std::cout << pow((integer)_p_,_d_) << std::endl;
			linbox_check((p ==  pow((integer)_p_,_d_)) );
			// std::cout << _p_ << std::endl;
			// typedef typename Field::Element Element ;
			typedef double Element;

			const BlasMatrix<Element > Zerooo(_r_,_s_); // add id ?
			std::vector<BlasMatrix<Element > > NN(_d_,Zerooo); // le noyal !!

			RightPadicNullspace((std::vector<BlasMatrix<Element > >&) NN,
					    (const BlasMatrix<Integer> &)_C_,
					    (const BlasMatrix<Integer> &)_B_,
					    _p_,
					    (int)_d_,
					    _meth_);

			BlasMatrix<Integer> K(_n_,_s_);
			N = Dadize(K,NN,_p_);
			Integer pn = Integer(1)-p;
			pn /= (Integer)(_p_-1);
			for (size_t i = 0; i < _s_; i++)
				N.setEntry( _r_+i , i , pn) ;


			// N.write(std::cout << "N_" << p << _d_ << ":= ",true) << ';' << std::endl;
			return N ;

		}

		// BlasMatrix<Element> initialize()
		// {
		// return BlasMatrix<Element>(m_,n_);
		// }

	};

	template<class Field> struct CRATemporaryVectorTrait< NS_mod_pd<Field>, typename Field::Element>
	{
		//typedef BlasMatrix<double>::pointer Type_t;
		typedef BlasMatrix<typename Field::Element> Type_t;
		// typedef FixedBlasMatrixDouble Type_t;
	};


	struct GenPowerPrime {
		Integer  _p_ ;
		long     _d_ ;
		typedef Integer Prime_Type ;
		RandomPrimeIterator _Pit_ ;
		GenPowerPrime( int _s , long d ) :
			_d_(d), _Pit_(RandomPrimeIterator(_s))
		{}
		inline GenPowerPrime &operator ++ ()
		{
			_p_ = pow(*(++_Pit_),_d_) ;
			return *this;
		}

		/** @brief get the random prime.
		 *  returns the actual prime.
		 */
		const Prime_Type &operator *  () const
		{
			return _p_;
		}

		/** @brief get the random prime.
		 *  returns the actual prime.
		 *  @warning a new prime is not generated.
		 */
		const Prime_Type & randomPrime() const
		{
			return _p_;
		}

		/** @brief Sets the seed.
		 *  Set the random seed to be \p ul.
		 *  @param ul the new seed.
		 */
		void static setSeed(unsigned long ul)
		{
			Integer::seeding(ul);
		}

	};

	template<class IMatrix>
	int RightCradicNullspace(IMatrix & N, const IMatrix & A, NullSpaceMethod   meth)
	{
		size_t m = A.rowdim();
		size_t n = A.coldim();
		size_t r = 0  ;    // rank
		size_t s ;         // ker dim, the result.

		typedef Modular<double> Field;     // dans lequel on réduit A
		typedef Field::Element  Element ;

		PID_integer ZZ ;
		RandomPrimeIter PrimeGen(19) ;
		Element p = (Element) PrimeGen.random();

		Field F(p);

		typedef Modular<Integer> Zpz ;

		MatrixDomain<PID_integer> ZMD(ZZ);

		BlasPermutation<size_t> P (n);
		BlasPermutation<size_t> Qt(m);

		{ /* LU mod p  */
			//!@todo reuse this info !!
			BlasMatrix<Element> Ap(m,n);
			MatrixHom::map(Ap,A,F);
			LQUPMatrix<Field> LUp(F,Ap,P,Qt);
			r = LUp.getrank();
		}

		s = n - r; // taille du noyal.
		// std::cout << "ker dim is (probably) " << s << std::endl;
		// d.resize(s);

		if (s == 0) { // A était inversible !
			// std::cout << "invertible !" << std::endl;
			N= BlasMatrix<Integer>(0,0);
			/* A is Probably invertible. second chance ? */
			// std::cout << "not checking again" << std::endl;
			return s ;
		}
		else if (r == 0) { // A est nulle ?
			// std::cout << "nulle !" << std::endl;
			bool nulle = ZMD.isZero(A) ;
			if (!nulle) {
				std::cout << "error : bad p !" << std::endl;
				return -1 ; // fail
			}
			N = BlasMatrix<Integer>(n,n);
			Integer one(1);
			for (size_t i = 0 ; i < n ; ++i){
				N.setEntry(i,i,one);
			}
			// std::cout << "ok, kernel is Id" << std::endl;
			return s ;
		}
		else { // A a un noyau non trivial !
			// std::cout << "on a un noyal !" << std::endl;

			/* pre-allocating it */

			/*  extracting C/B */
			BlasMatrix<Integer> * C = NULL ;
			BlasMatrix<Integer> * B = NULL ;

			P.Compress();
			Qt.Compress();

			if (!Qt.isIdentity()) {
				Qt.Invert();
				Qt.Compress();
			}
			if (P.isIdentity() && Qt.isIdentity()) {
				C = (new BlasMatrix<Integer>(A,0,0,r,r));
				B = (new BlasMatrix<Integer>(A,0,r,r,s));
			}
			else {
				C = (new BlasMatrix<Integer>(r,r)) ;
				B = (new BlasMatrix<Integer>(r,s)) ;
			}

			StackMat<BlasMatrix<Integer> >  QAP(*C,*B) ; // A = (C|B\?|?)
			if (!P.isIdentity() || !Qt.isIdentity()) {
				extractMatrix(QAP,A,Qt,P) ;
			}
			size_t lA ;
			MaxElementSize(lA,A) ;

			/*  Nullspace mod pi^d */
			unsigned int d = 20 ;
			size_t PrimeSize  = 10;
			double UpperBound = (_LB_LOG2 * (double) lA + log(sqrt(double(n))))*(double) n ;
			GenPowerPrime genprime( PrimeSize, d );
			// RandomPrimeIterator genprime( PrimeSize );

			ChineseRemainder< FullMultipBlasMatCRA< Zpz > > cra( std::pair<size_t,double>(n*s, UpperBound) );

			NS_mod_pd<Zpz> Iteration(*C,*B,d,meth);

			BlasMatrix<Integer> NN(n,s);
			cra(NN, Iteration, genprime);
			// BlasMatrix<Integer> NNN(n,s);
			N = NN ;
			//!@todo submatrices !
			// for (size_t i = 0 ; i < r ; ++i)
				// for (size_t j = 0 ; j < s ; ++j)
					// N.setEntry(i,j,NN.getEntry(i,j));
			/*  aggregating parts  */

			// write(std::cout << "Z := ",NN,p,true) << ";" << std::endl;

			// Integer pn = Integer(1)-pow((Integer)p,(long)d);
			// pn /= (Integer)(p-1);
			// for (size_t i = 0; i < s; i++)
				// N.setEntry( r+i , i , pn) ;
			BlasMatrixDomain<PID_integer> ZZDom(ZZ);
			TransposedBlasMatrix<BlasPermutation<size_t > > Pt(P);
			ZZDom.mulin_right (Pt,N);

			// N.write(std::cout << "N := ",true) << ';' << std::endl;


			/*-------------------*/
		}


		return s ;

	}

#if 0
	template<class Randiter>
	BlasMatrix<Integer> &
	random_lu_rank(const PID_integer & F,
		       const Randiter & R,
		       BlasMatrix<Integer> &A,
		       int & rank)
	{

		typedef Modular<double> Field ;
		typedef Field::Element Element ;
		/*- CRA -*/
		// ChineseRemainder< FullMultipFixedCRA< Field > > cra( std::pair<size_t,double>(A.rowdim()*A.coldim(), UpperBound) );


		// typename BlasMatrix<Integer>::RawIterator A_it = A.rawBegin() ;


		return A ;
	}
#endif
#endif
}

#undef _LB_MAX_NS_TRY
#undef _LB_LOG2

#endif //__LINBOX_integer_nullspace_INL


